Spatial Accessibility

Taken from: Spatial Methods in Community Research by Professor Noli Brazil https://crd230.github.io/lab5a.html#

install.packages("matrixStats")
install.packages("SpatialAcc")

** SDAD Note: Before some packages will load correctly we need to make sure V8 is installed correctly. To do so, run the following: Sys.setenv(DOWNLOAD_STATIC_LIBV8 = 1) install.packages(“V8”)

library(sf)

Linking to GEOS 3.7.2, GDAL 3.0.0, PROJ 6.1.1

library(tidyverse)

── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.2     ✓ dplyr   1.0.7
✓ tidyr   1.0.3     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.0
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(tmap)
library(tidycensus)
library(tigris)

To enable 
caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.

library(rmapshaper)

Registered S3 method overwritten by 'geojsonlint':
  method         from 
  print.location dplyr

library(matrixStats)

Attaching package: ‘matrixStats’

The following object is masked from ‘package:dplyr’:

    count

library(SpatialAcc)

#load in your Census API Key
census_api_key(Sys.getenv("CENSUS_API_KEY"))

Read in census tract data

We will need to bring in census tract polygon features and racial composition data using the Census API and keep tracts within Los Angeles city boundaries. The code for accomplishing these tasks is below. We won’t go through each line of code in detail because we’ve covered all of these operations and functions in prior labs. I’ve embedded comments within the code that briefly explain what each chunk is doing, but go back to prior guides (or RDS/GWR) if you need further help. Note the use of output = wide in get_acs() to start out with a wide rather than a long dataset.

#Clip tracts using LA boundary
la.city.tracts <- ms_clip(target = ca.tracts, clip = la.city, remove_slivers = TRUE)

Error in ms_clip(target = ca.tracts, clip = la.city, remove_slivers = TRUE) : 
  could not find function "ms_clip"

Read in the bank data

Let’s next bring in a shapefile of full-service banks and credit unions in Los Angeles as of 2020. I zipped up the file and uploaded it onto Github. Use the following code to download and unzip the file.

download.file(url = "https://raw.githubusercontent.com/crd230/data/master/los_angeles_banks.zip", destfile = "los_angeles_banks.zip")

trying URL 'https://raw.githubusercontent.com/crd230/data/master/los_angeles_banks.zip'
Content type 'application/zip' length 27829 bytes (27 KB)
==================================================
downloaded 27 KB

unzip(zipfile = "los_angeles_banks.zip")

Don’t worry if you don’t understand these commands - they are more for you to simply copy and paste so that you can download files that I zipped up and uploaded onto GitHub. You can look at the help documentation for each function if you are curious. Bring the data into R using st_read().

banks <-st_read("los_angeles_banks.shp")

Reading layer `los_angeles_banks' from data source `/home/ads7fg/git/VDH/src/los_angeles_banks.shp' using driver `ESRI Shapefile'
Simple feature collection with 616 features and 7 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 348646.4 ymin: 3733033 xmax: 392442.4 ymax: 3797503
CRS:            unknown

The data are already in a spatial format, but if they came in street addresses or latitudes and longitudes, you would follow the methods described in Lab 4a to convert them to a spatial format.

Bank branch data come from the Federal Deposit Insurance Corporation. A bank branch is a physical location of a banking corporation, such as Chase, Bank of America or Wells Fargo. These buildings are technically referred to as “brick-and-mortar” branches, and they provide face-to-face service for customers of a bank. Credit union branch data come from the National Credit Union Administration. A credit union is a nonprofit financial institution that’s owned by the people who use its financial products. Credit union members can access the same kinds of products and services as offered by a traditional bank, such as credit cards, checking and savings accounts and loans. Members elect a board of directors to manage the credit union to ensure that their best interests are represented.

Take a look at the data

glimpse(banks)

Rows: 616
Columns: 8
$ name     <chr> "MUFG Union Bank, National Association", "MUFG Union Bank, National Association", "MUFG Union Bank, National Association", "MUFG Union Bank, National Association", "MUFG Union Bank, National Asso…
$ Address  <chr> "120 South San Pedro Street", "3060 Crenshaw Boulevard", "14360 Roscoe Boulevard", "5855 Topanga Canyon Blvd", "16633 Ventura Boulevard", "101 South San Pedro Street", "3810 Crenshaw Boulevard", …
$ City     <chr> "Los Angeles", "Los Angeles", "Panorama City", "Woodland Hills", "Encino", "Los Angeles", "Los Angeles", "Los Angeles", "Los Angeles", "Pacoima", "North Hollywood", "Los Angeles", "Los Angeles", …
$ State    <chr> "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA…
$ zipcode  <chr> "90012", "90016", "91402", "91367", "91436", "90012", "90008", "90064", "90010", "91331", "91606", "90024", "90048", "90067", "91436", "91436", "91364", "90026", "91352", "91316", "90005", "91367…
$ deposits <dbl> 279.954, 63.514, 112.858, 289.353, 234.795, 149.829, 51.791, 154.548, 347.492, 62.000, 91.000, 245.885, 132.697, 709.914, 290.002, 96.391, 126.372, 23.010, 83.889, 61.777, 212.616, 90.722, 300.80…
$ creditu  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ geometry <POINT [m]> POINT (385459.1 3768403), POINT (376750.7 3765979), POINT (366769.1 3787663), POINT (352003.6 3782627), POINT (362131.5 3780698), POINT (385458.6 3768433), POINT (376754.3 3764717), POINT (…

Make sure the tract and bank Coordinate Reference Systems are the same and in UTM.

st_crs(la.city.tracts.utm)$proj4string

[1] "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"

st_crs(banks)$proj4string

[1] "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"

Let’s map the banks on top of the city’s tracts using our pal tm_shape()

Points in polygon

In its most basic form, spatial accessibility factors in the supply of an amenity, the demand for that supply, and the distance facilitating or impeding the linkage between supply and demand. A simple approach to measuring spatial access is to count the number of amenities located within neighborhood boundaries. We already covered how to do this in Lab 4a. To get the count of banks by census tract, we use st_join() to perform a spatial join of the banks to the tracts and then group_by() and sumarize() to execute the count.

la.city.tracts.banks

Simple feature collection with 234 features and 2 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 348646.4 ymin: 3733033 xmax: 392442.4 ymax: 3797503
CRS:            unknown

We’ll join la.city.tracts.banks back to the original tract file. We have to first drop la.city.tracts.banks’s geometry because you cannot join two sf objects together using left_join().

Now, use left_join() and assign a 0 to tracts with no banks using the replace_na() function within mutate().

The number of banks represents the supply. Distance is captured by banks residing within tract boundaries. The demand can take on many forms, but one common proxy is the resident population size. Let’s calculate the number of banks per resident population, multiplying the ratio by 10,000 to get more meaningful accessibility ratio sizes.

What does the distribution of this variable look like?

summary(la.city.tracts.utm$banksperpop)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   0.00    0.00    0.00    1.78    0.00  105.26       4 

Let’s create a choropleth map of banks per 10,000 population.

Point proximity buffers

The points-in-polygon approach treats neighborhood boundaries as impermeable as if banks right outside boundary lines are unreachable to residents. Another approach to measuring access is to capture distance from the nearest amenities. There are several ways to account for distance. One method involves counting the number of amenities within r miles of the neighborhood. The number of amenities within r refers to the cumulative opportunity, or, in other words, to the availability provided by the immediate surroundings.

The process involves generating a circle with a radius r (i.e. a circular buffer of size r) around existing geographic features. You connect other features to the points based on whether they fall inside or outside the boundary of the buffer. The buffer defines this filtering window as that location’s catchment area. Accessibility is measured as the ratio of amenities to residents (or supply to demand) located within the buffer.

The first task is to create buffers of radius r around each neighborhood. We draw a buffer around the neighborhood centroid, which is the point on the middle of a tract. To generate centroids, use the sf function st_centroid()

tract.centroids <- st_centroid(la.city.tracts.utm)

st_centroid assumes attributes are constant over geometries of x

Plot the centroids on top of the tracts so we know what we produced

We then create circular buffers of distance r around each tract’s centroid. The size of the buffer radius depends on the city you are examining and the context of your question. To start, let’s use 1-mile or 1609-meter buffers.

We use the sf function st_buffer() to create buffers. The required arguments are your sf object and the distance dist =. Remember that the units for UTM is meters, so specifying dist = 1609 means 1609 meters.

tract.buff  <-st_buffer(tract.centroids, dist = 1609)
tract.buff

Simple feature collection with 1001 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 345766.8 ymin: 3729537 xmax: 394269.8 ymax: 3800782
CRS:            +proj=utm +zone=11N +datum=NAD83 +ellps=GRS80
First 10 features:
         GEOID tpop  pnhwhite    pnhasn    pnhblk     phisp medinc banks                       geometry banksperpop
1  06037199800 5828 0.0130405 0.3637612 0.0001716 0.6230268  37755     0 POLYGON ((389925.6 3770651,...    0.000000
2  06037291300 3037 0.1471847 0.4333224 0.0951597 0.2462957  81281     0 POLYGON ((381974.1 3747942,...    0.000000
3  06037292000 6567 0.0753769 0.1770976 0.1067459 0.5541343  42135     0 POLYGON ((381348.7 3746245,...    0.000000
4  06037209510 3076 0.0325098 0.2334200 0.0159298 0.6849805  30722     0 POLYGON ((383292.7 3768770,...    0.000000
5  06037267501 5993 0.4768897 0.1558485 0.0258635 0.2848323  86086     4 POLYGON ((367206.7 3767693,...    6.674454
6  06037294520 3710 0.0353100 0.0064690 0.0215633 0.9342318  50732     0 POLYGON ((384418.5 3739241,...    0.000000
7  06037269904 3927 0.5059842 0.1807996 0.1245225 0.1451490  70682     0 POLYGON ((371264.4 3766137,...    0.000000
8  06037183820 3445 0.0998549 0.0928882 0.0272859 0.7582003  46477     0 POLYGON ((391198.1 3774332,...    0.000000
9  06037197300 4052 0.4437315 0.1591807 0.0338105 0.3284798  83843     0 POLYGON ((386086.8 3771868,...    0.000000
10 06037204920 2751 0.0170847 0.0000000 0.0083606 0.9727372  39904     0 POLYGON ((391042.1 3764745,...    0.000000

You’ll see that tract.buff is a polygon object like a census tract. To be clear regarding what a buffer is, let’s extract one tract, GEOID = 06037262501 , and its 1-mile buffer.

ex1 <- filter(tract.buff, GEOID == "06037262501")
ex2 <- filter(la.city.tracts.utm, GEOID == "06037262501")

And let’s map it onto tracts with the banks.

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(la.city.tracts.utm) +
          tm_polygons() +
        tm_shape(tract.centroids) +
          tm_dots(size = 0.01) +
          tm_shape(ex1) +
          tm_borders(col="red") +  
        tm_shape(ex2) +
          tm_dots(col = "red") 

The tract we are mapping is located in the West Los Angeles area along the coast. If you zoom in, you’ll see that the centroid is right in the middle of the tract and the buffer is a one-mile circle surrounding the tract. The buffer contains not just tract 06037262501, but touches other nearby tracts.

We want to count the number of banks located within the buffers. Remember that tract.buff is a polygon. Counting the number of banks inside a buffer is a points-in-polygons operation, which we already did above. We’ll calculate the number of banks per 10,000 resident population.

buff.banks <- banks  %>% 
        st_join(tract.buff) %>%
        group_by(GEOID) %>% 
        summarize(banks1m = n()) 

buff.banks <- st_drop_geometry(buff.banks)

la.city.tracts.utm <- la.city.tracts.utm %>%
                      left_join(buff.banks, by = "GEOID") %>%
                      mutate(banks1m = replace_na(banks1m, 0)) %>%
                      mutate(banksbuff1m = (banks1m/tpop)*10000)

We can then use summarize() to find the mean number of banks that occur within 1 mile of a tract.

There are on average 6.9 banks within 1 mile of tracts in Los Angeles. This is likely not representative given that there are a small number of tracts with a lot of banks within a 1-mile radius, as demonstrated by the histogram below (the tract with 70+ banks is in the financial district). It’s common practice to test the sensitivity of your results by choosing different buffer sizes.

Let’s create a choropleth map of banks per population within a 1-mile radius.

tmap_mode("plot")

tmap mode set to plotting

tm_shape(la.city.tracts.utm, unit = "mi") +
  tm_polygons(col = "banksbuff1m", style = "jenks",palette = "Reds", 
              border.alpha = 0, title = "Banks per\n10k population") +
  tm_scale_bar(position = c("left", "bottom")) +
    tm_layout(main.title = "Bank spatial accessibility in Los Angeles City",
            main.title.size = 0.95, frame = FALSE,
            legend.outside = TRUE, legend.outside.position = "right")

Distance to nearest bank

We limit ourselves a bit by capturing banks just within a radius distance. Why not calculate the distance to all banks and summarize that distribution? Distance is typically measured in Euclidean distance. Euclidean distance, also referred to as straight-line distance or “as the crow flies,” is the distance between two points connected by a straight line on a flat surface. To calculate the distance (in meters) from each tract’s centroid to each bank we use the sf function st_distance().

For the object bank.dist, the rows represent the tracts and the columns are the banks. You can check this by comparing dimensions.

#number of tracts is 999
dim(tract.centroids)

[1] 1001   10

#number of banks is 622
dim(banks)

[1] 616   8

#999 by 622
dim(bank.dist)

[1] 1001  616

What we have with bank.dist is a distance matrix. That is, we have each tracts’s distance to each bank. Note that st_distance() will calculate distance using the Coordinate Reference System’s units, in this case meters. A common spatial accessibility measure is the distance to the closest amenity (e.g. nearest bank). It allows one to evaluate the immediate proximity to the bank. We can use the function rowMins() in the package matrixStats to calculate the shortest distance.

The function rowMins() does exactly what you think it would do - get the minimum value across columns for each row. For bank.dist, this means we get the minimum distance to a bank for each neighborhood. Note that the order of tracts in bank.dist is the same as in la.city.tracts.utm. Therefore, we can run rowMins() within mutate() to save the resulting value in our main data set la.city.tracts.utm.

We can use summarize() to find the mean minimum distance to banks.

The closest bank to a neighborhood in Los Angeles is on average 1,136 meters, or about 3/4 a mile.

Let’s create a choropleth map of nearest bank distance.

Rather than calculating the distance to the nearest amenity, another approach is to measure the mean distance to the n closest services. This approach conceptualizes accessibility as the average cost to reach a diversity of options rather than a single option.

Floating Catchment Area

The Floating Catchment Area (FCA) method extends the buffer approach by incorporating the demand of other nearby tracts for an amenity. Similar to the basic buffer approach, FCA draws a circle (or buffer) of radius d around a tract’s centroid and defines this filtering window as that location’s catchment area. In contrast to the simple buffer method, FCA accounts for not just the number of amenities that fall inside the buffer, but also the centroids of other nearby tracts. The method accounts for competition for services among residents in other neighborhoods (competing demand).

We’ve already calculated the FCA numerator in the buffer analysis (banks1m). We apply the buffer analysis to count tract centroids tract.centroids rather than banks within a 1 mile buffer of each tract tract.buff. We sum up the populations of tracts whose centroids fall in each buffer.

buff.tracts <- tract.centroids  %>% 
        select(tpop) %>% #we don't need other variables
        st_join(tract.buff) %>%
        group_by(GEOID) %>% 
        summarize(buffpop = sum(tpop.x)) %>%
        ungroup()

buff.tracts <- st_drop_geometry(buff.tracts)

The variable buffpop in buff.tracts reflects the total population of tracts whose centroids fall into a tract’s 1-mile buffer. We join this variable back into our main tibble.

We next create the variable fca by dividing the number of banks in the buffer’s tract by the tract’s population and the population of nearby tracts whose centroids fall inside the buffer. We multiply by 10,000 to make it a per 10,000 population rate.

Map it. Higher values indicates greater spatial accessibility.

Two-step Floating Catchment Area

The basic Floating Catchment Area method calculates access strictly from the demand (resident population) perspective. The method does not take into account the service area of the supplier (bank). The two-step Floating Catchment Area (2SFCA) method remedies this issue by calculating catchment areas for both supplier and consumer.

To calculate 2SFCA ratios, we use the function ac() which is part of the SpatialAcc package. There are five inputs you need to supply ac(). First, a matrix of distances between tracts to banks. We already did that above when we constructed bank.dist, but the resulting object is not a matrix, but has class “units”

class(bank.dist)

[1] "units"

Let’s instead use the function distance() which is part of the SpatialAcc package. In order to use this function, you will need to first get the coordinates of the tract centroids and banks. Use the sf function st_coordinates().

Then plug these into the function distance(), specifying type = “euclidean” to calculate euclidean distances.

class(dist.matrix)

[1] "matrix" "array" 

In addition to the distance matrix, we also need to specify measures of demand p and supply n in the function ac(). The measures we’ve used so far proxy the demand for banks as total resident population. What about supply? In the basic Floating Catchment Area method, we treated every bank equally. In other words, we assumed that capacity and quality are equal across banks. This is likely not true. One measure of bank capacity (less about quality) is the number of its deposits in dollars. This is measured by the variable deposits in the file banks, which measures the amount of deposits in millions. Plug these demand and supply values into ac() to get 2SFCA measures of accessibility.

We also had to specify the threshold distance that defines the catchment area (like the buffer radius r). Here, we use 1609 meters, or one mile. The resulting object is a vector of 2SFCA values. Save the results in la.city.tracts.utm.

Then mapificate

Higher values of TSFCA equals greater spatial access.